###免疫细胞热图

#引用包
library(limma)
library(pheatmap)

clusterFile="rs_tcga.txt"      #分型的结果文件
immFile="infiltration_estimation_for_tcga.csv"     #免疫细胞浸润结果文件
setwd("E:\\1.Ovarian_Treg\\1.picture\\11.IMMUNE")      #设置工作目录

#读取分型结果文件
cluster=read.table(clusterFile, header=T, sep="\t", check.names=F, row.names=1)
cluster$risk=ifelse(cluster$riskScore>median(cluster$riskScore),"high","low")
#读取免疫细胞浸润文件
immune=read.csv(immFile, header=T, sep=",", check.names=F, row.names=1)
immune=as.matrix(immune)
rownames(immune)=gsub("(.*?)\\-(.*?)\\-(.*?)\\-(.*)", "\\1\\-\\2\\-\\3", rownames(immune))
immune=avereps(immune)

#病人分型和免疫细胞合并
sameSample=intersect(row.names(cluster), row.names(immune))
cluster=cluster[sameSample, "risk", drop=F]
immune=immune[sameSample, , drop=F]
data=cbind(cluster, immune)

#分型的免疫细胞差异分析
outTab=data.frame()
sigCell=c("risk")
for(i in colnames(data)[2:ncol(data)]){
  if(sd(data[,i])<0.001){next}
  if(length(levels(factor(data[,"risk"])))>2){
    test=kruskal.test(data[,i] ~ data[,"risk"])
  }else{
    test=wilcox.test(data[,i] ~ data[,"risk"])
  }
  pvalue=test$p.value
  if(pvalue<0.05){
    outTab=rbind(outTab,cbind(immune=i, pvalue))
    sigCell=c(sigCell, i)
  }
}
write.table(file="immuneCor.txt", outTab, sep="\t", quote=F, row.names=F)

#热图数据
data=data[,sigCell]
#分型的注释
data=data[order(data[,"risk"]),]
annCol=data[,1,drop=F]
annCol[,"risk"]=factor(annCol[,"risk"], unique(annCol[,"risk"]))
data=t(data[,(2:ncol(data))])
#软件类型的注释
annRow=sapply(strsplit(rownames(data),"_"), '[', 2)
annRow=as.data.frame(annRow)
row.names(annRow)=row.names(data)
colnames(annRow)=c("Methods")
annRow[,"Methods"]=factor(annRow[,"Methods"], unique(annRow[,"Methods"]))
gapCol=as.vector(cumsum(table(annCol[,"risk"])))
gapRow=as.vector(cumsum(table(annRow[,"Methods"])))

#定义热图注释的颜色
bioCol=ggsci::pal_npg("nrc")(10)
bioCol=bioCol[1:length(unique(annCol[,"risk"]))]
risk=bioCol
names(risk)=levels(factor(annCol[,"risk"]))

##
bioCol=ggsci::pal_npg("nrc")(10)
bioCol=bioCol[3:10][1:length(unique(annRow[,"Methods"]))]
Methods=bioCol
names(Methods)=levels(factor(annRow[,"Methods"]))
ann_colors=list(risk=risk,Methods=Methods)

#热图可视化



p1=pheatmap(data,
            annotation=annCol,
            annotation_row=annRow,
            annotation_colors = ann_colors,
            color = colorRampPalette(c(rep("#2F62A0",3), "white", rep("#920000",3)))(100),
            cluster_cols =F,
            cluster_rows =F,
            gaps_row=gapRow,
            gaps_col=gapCol,
            scale="row",
            show_colnames=F,
            show_rownames=T,
            fontsize=6,
            fontsize_row=5,
            fontsize_col=6)


